Salamander eye scaling and investment

Data

First we import our tidied salamander morphology & trait dataset.

# Import tidied salamander eye size and trait data
salamanders <- data.frame(read.csv("../Data/Tidy/salamanders_tidy.csv", header=TRUE, na.strings=c("", "NA", " ","<1")))

# Quick look at data structure
str(salamanders)

## 'data.frame':    424 obs. of  39 variables:
##  $ X                 : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Genus_Species     : chr  "Ambystoma_amblycephalum" "Ambystoma_amblycephalum" "Ambystoma_annulatum" "Ambystoma_annulatum" ...
##  $ Order             : chr  "Caudata" "Caudata" "Caudata" "Caudata" ...
##  $ Suborder          : chr  "Salamandroidea" "Salamandroidea" "Salamandroidea" "Salamandroidea" ...
##  $ Family            : chr  "Ambystomatidae" "Ambystomatidae" "Ambystomatidae" "Ambystomatidae" ...
##  $ Subfamily         : chr  NA NA NA NA ...
##  $ Genus             : chr  "Ambystoma" "Ambystoma" "Ambystoma" "Ambystoma" ...
##  $ Species           : chr  "amblycephalum" "amblycephalum" "annulatum" "annulatum" ...
##  $ SVL_mm            : int  73 60 80 80 88 60 56 102 103 115 ...
##  $ Mass_g            : num  11.58 9.34 8.97 8.39 12.12 ...
##  $ rootmass          : num  2.26 2.11 2.08 2.03 2.3 ...
##  $ ED_right_mm       : num  4.11 3.82 3.1 3.45 3.4 2.98 2.24 5.25 4.77 4.95 ...
##  $ ED_left_mm        : num  3.96 3.63 3.17 3.51 3.6 2.58 2.17 5.2 4.83 4.87 ...
##  $ eyemean           : num  4.04 3.73 3.13 3.48 3.5 ...
##  $ CD_right_mm       : num  2.93 2.42 2.06 2.34 2.2 1.7 1.52 3.37 3.56 3.62 ...
##  $ CD_left_mm        : num  2.94 2.11 1.99 2.44 2.43 1.7 1.75 3.32 3.33 4.11 ...
##  $ cormean           : num  2.94 2.27 2.02 2.39 2.31 ...
##  $ Gill_Presence     : chr  "yes" "yes" "no" "no" ...
##  $ Aquatic           : chr  "yes" "yes" "no" "no" ...
##  $ Semiaquatic       : chr  "no" "no" "no" "no" ...
##  $ Terrestrial       : chr  "no" "no" "no" "no" ...
##  $ Scansorial        : chr  "no" "no" "no" "no" ...
##  $ Fossorial         : chr  "no" "no" "yes" "yes" ...
##  $ Subfossorial      : chr  "no" "no" "no" "no" ...
##  $ Cave_dwelling     : chr  "no" "no" "no" "no" ...
##  $ Nocturnal         : chr  NA NA "yes" "yes" ...
##  $ Diurnal           : chr  NA NA "no" "no" ...
##  $ Paedomorphic      : chr  "yes" "yes" "no" "no" ...
##  $ Metamorphosizing  : chr  "no" "no" "yes" "yes" ...
##  $ Direct_Development: chr  "no" "no" "no" "no" ...
##  $ Free_living       : chr  "yes" "yes" "yes" "yes" ...
##  $ No_free_living    : chr  "no" "no" "no" "no" ...
##  $ Lentic_water      : chr  "yes" "yes" "yes" "yes" ...
##  $ Lotic_water       : chr  "no" "no" "no" "no" ...
##  $ Both_water        : chr  "no" "no" "no" "no" ...
##  $ No_larva          : chr  "no" "no" "no" "no" ...
##  $ Females_larger    : chr  "yes" "yes" "yes" "yes" ...
##  $ Males_larger      : chr  "no" "no" "no" "no" ...
##  $ No_SSD            : chr  "no" "no" "no" "no" ...

Next, we import the amphibian tree from Jetz and Pyron (2019) that has been pruned to match the salamander species in this dataset.

#Import pruned phylogeny
caudatatree <- read.nexus("../Data/Tidy/caudata-tree.nex")

#Plot tree
plot.phylo(caudatatree, show.tip.label=TRUE, cex = 0.7)

Data exploration

First, we can take a preliminary look at all of our data so that we can see how it looks and whether we notice any obvious outliers that we might want to go back and check out. Here we plot out our measurements from the right side vs. the left side of each specimen so that we can make sure that we have consistent measurements throughout the dataset. We use the package plotly to create interactive plots - hover over a point to see the name of the species measured.

Right eyes vs. left eyes

We expect that measurements of left and right eyes and corneas should be similar within specimens. We can plot out our data to test for this and make sure we don’t see any glaring issues. We can also fit a regression, which we would expect to have a slope near 1 and an intercept near 0 if left and right eye measurments are similar.

####right eye diameter vs. left eye diameter ####

#linear model (ordinary least squares model)
RLeye.lm <- lm(ED_right_mm ~ ED_left_mm, data = salamanders)
summary(RLeye.lm)

## 
## Call:
## lm(formula = ED_right_mm ~ ED_left_mm, data = salamanders)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53420 -0.13491 -0.01135  0.12950  0.67196 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.015122   0.035243   0.429    0.668    
## ED_left_mm  1.010149   0.009162 110.253   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2045 on 383 degrees of freedom
##   (39 observations deleted due to missingness)
## Multiple R-squared:  0.9695, Adjusted R-squared:  0.9694 
## F-statistic: 1.216e+04 on 1 and 383 DF,  p-value: < 2.2e-16

#plot
RLeye.plot <- ggplot(salamanders, aes(x = ED_left_mm, y = ED_right_mm, text = Genus_Species)) +
 geom_point(alpha = 0.9) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(RLeye.lm)[[2]], intercept = coef(RLeye.lm)[[1]]) #plots OLS fit 


#interactive plot
ggplotly(RLeye.plot)

As we would expect, measurements of the left and right eye are highly correlated, with a slope close to 1 and an intercept close to 0. We can do the same test for the corneal measurements.

####right cornea diameter vs. left cornea diameter ####

#linear model (ordinary least squares model)
RLcornea.lm <- lm(CD_right_mm ~ CD_left_mm, data = salamanders)
summary(RLcornea.lm)

## 
## Call:
## lm(formula = CD_right_mm ~ CD_left_mm, data = salamanders)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71444 -0.11655 -0.00936  0.13772  0.78095 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.10579    0.03248   3.257  0.00122 ** 
## CD_left_mm   0.97409    0.01177  82.742  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2133 on 422 degrees of freedom
## Multiple R-squared:  0.9419, Adjusted R-squared:  0.9418 
## F-statistic:  6846 on 1 and 422 DF,  p-value: < 2.2e-16

#plot
RLcor.plot <- ggplot(salamanders, aes(x = CD_left_mm, y = CD_right_mm, text = Genus_Species)) +
 geom_point(alpha = 0.9) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(RLcornea.lm)[[2]], intercept = coef(RLcornea.lm)[[1]]) #plots OLS fit 


#interactive plot
ggplotly(RLcor.plot)

Again, we find that left and right corneal measurements are highly correlated, with a slope close to 1 and an intercept close to 0. Together, these results indicate that measurements were consistent between the left and right eyes of a specimen, which is good! For looking at interspecific allometry, we will calculate species means from these data so that we can use phylogenetic comparative methods.

Interspecific allometry

Here, we will use phylogenetic least-squares (PGLS) regressions to see how eyes are scaling with body size across different species of salamanders with different ecologies. PGLS analyses are similar to normal linear regressions, but also account for the non-independence of data collected from close relatives due to shared evolutionary history.

In these analyses, we will log eye size and body size variables, as this is standard for allometric relationships.

Species means for morphological comparisons

First, we will calculate the mean measurements across specimens for each species. We need to make sure in doing this that we omit rows with missing data for a variable of interest. In looking at the data, we can see that we have many more measurements for cornea diameter than for eye diameter

#find species with missing eye measurements
eye.missing <- salamanders %>%
  filter(is.na(eyemean)) %>%
  droplevels() %>%
  group_by(Family,Genus_Species) %>%
  summarise(n = n())

#create RMarkdown table of species missing eye data
kable(eye.missing, caption = "Specimens without eye size data") %>%
  kable_styling(full_width = F) %>%
  collapse_rows(columns = 1, valign = "top") %>%
  scroll_box(height = "500px")

Specimens without eye size data

Family	Genus_Species	n
Amphiumidae	Amphiuma_means	3
Amphiumidae	Amphiuma_tridactylum	3
Cryptobranchidae	Andrias_davidianus	2
Cryptobranchidae	Andrias_japonicus	3
Cryptobranchidae	Cryptobranchus_alleganiensis	4
Plethodontidae	Eurycea_nana	3
Plethodontidae	Eurycea_rathbuni	3
Proteidae	Necturus_alabamensis	1
Proteidae	Necturus_beyeri	2
Proteidae	Necturus_lewisi	3
Proteidae	Necturus_maculosus	3
Salamandridae	Pseudobranchus_striatus	3
Sirenidae	Siren_intermedia	3
Sirenidae	Siren_lacertina	3

#find species missing cornea measurements
cor.missing <- salamanders %>%
  filter(is.na(cormean)) %>%
  droplevels() %>%
  group_by(Family,Genus_Species) %>%
  summarise(n = n())

#species missing cornea measurements
print(cor.missing)

## # A tibble: 0 × 3
## # Groups:   Family [0]
## # … with 3 variables: Family <chr>, Genus_Species <chr>, n <int>
## # ℹ Use `colnames()` to see all variable names

We then find species means for various morphological parameters. There are some specimens that lack data for eye size but have data for svl, mass, and corneal diameter (these measurements are present for all specimens). However, when eye size measurements are missing, they are consistently missing for all specimens measured within a species. This means we don’t have to worry about missing data for one specimen corrupting the mean across several specimens for eye size (it is either present for all specimens in a species or absent for all specimens in a species). So we can make just one dataframe of species averages, and some species will just lack eye size data.

# species means for cornea data
av.sal <- salamanders %>% 
  mutate_if(is.character, as.factor) %>%
  group_by(Genus_Species) %>%
  summarise(eye_av = mean(eyemean), 
            cor_av = mean(cormean),
            svl_av = mean(SVL_mm), 
            rootmass_av = mean(rootmass), 
            n = n())

## Merge data means with other column info, keeping only trait columns
av.sal <- merge(av.sal, salamanders[match(unique(salamanders$Genus_Species), salamanders$Genus_Species), ], by="Genus_Species", all.x = TRUE, all.y = FALSE) %>% 
  select(Genus_Species, eye_av, cor_av, svl_av, rootmass_av, n, Order, Suborder, Family, Subfamily, Genus, Species, Aquatic, Semiaquatic, Terrestrial, Scansorial, Fossorial, Subfossorial, Cave_dwelling, Nocturnal, Diurnal, Paedomorphic, Metamorphosizing, Direct_Development, Free_living, No_free_living, Lentic_water, Lotic_water, Both_water, No_larva, Females_larger, Males_larger, No_SSD) 

Phylogenetic generalised least squares regressions

In order to fit a PGLS regression in caper, we first need to make a comparative data object that includes our dataset and our phylogeny. Note that the gilled Pleurodeles waltl will be dropped from this, as there can only be one row for each species.

#Make row names of dataset the species names (so it will match phylogeny tips)
rownames(av.sal) <- av.sal$Genus_Species

#check that names match in dataframe and tree
name.check(phy = caudatatree, data = av.sal)

## $tree_not_data
## character(0)
## 
## $data_not_tree
## [1] "Pleurodeles_waltl_gilled"

#use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
sal.comp <- comparative.data(phy = caudatatree, data = av.sal, 
                            names.col = Genus_Species, vcv = TRUE, 
                            na.omit = FALSE, warn.dropped = TRUE)

#check for dropped tips or dropped species
sal.comp$dropped$tips #phylogeny

## character(0)

sal.comp$dropped$unmatched.rows #dataset

## [1] "Pleurodeles_waltl_gilled"

PGLS of log eye diameter vs. log cube root of mass

#PGLS eye diameter vs. cube root of mass
pgls_edmass <- pgls(log10(eye_av) ~ log10(rootmass_av), 
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edmass)

par(mfrow = c(1, 1))

#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.edmass <- pgls.profile(pgls_edmass, "lambda")
plot(lambda.edmass)

#print model output 
summary(pgls_edmass)

## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(rootmass_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0180704 -0.0036954 -0.0000669  0.0035026  0.0228160 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.845
##    lower bound : 0.000, p = 3.0531e-11
##    upper bound : 1.000, p = 2.8592e-06
##    95.0% CI   : (0.658, 0.942)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)        0.398691   0.042603  9.3584  2.22e-16 ***
## log10(rootmass_av) 0.773017   0.041458 18.6458 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.006815 on 136 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared: 0.7188,  Adjusted R-squared: 0.7167 
## F-statistic: 347.7 on 1 and 136 DF,  p-value: < 2.2e-16

#plot
edmass.plot <- ggplot(data = av.sal, aes(x = log10(rootmass_av), y = log10(eye_av), text = Genus_Species)) +
 geom_point(alpha = 0.7) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_edmass)[[2]], intercept = coef(pgls_edmass)[[1]]) #plots OLS fit 


#interactive plot
ggplotly(edmass.plot)

PGLS of log eye diameter vs. log snout-vent length

#PGLS eye diameter vs. SVL
pgls_edsvl <- pgls(log10(eye_av) ~ log10(svl_av), 
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edsvl)

par(mfrow = c(1, 1))

#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.edsvl <- pgls.profile(pgls_edsvl, "lambda")
plot(lambda.edsvl)

#print model output 
summary(pgls_edsvl)

## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(svl_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0173372 -0.0021347  0.0006661  0.0041418  0.0154655 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.749
##    lower bound : 0.000, p = 3.9968e-14
##    upper bound : 1.000, p = 1.1559e-08
##    95.0% CI   : (0.506, 0.898)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   -0.839890   0.079508 -10.564 < 2.2e-16 ***
## log10(svl_av)  0.782916   0.039652  19.745 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005834 on 136 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared: 0.7414,  Adjusted R-squared: 0.7395 
## F-statistic: 389.9 on 1 and 136 DF,  p-value: < 2.2e-16

#plot
edsvl.plot <- ggplot(data = av.sal, aes(x = log10(svl_av), y = log10(eye_av), text = Genus_Species)) +
 geom_point(alpha = 0.7) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_edsvl)[[2]], intercept = coef(pgls_edsvl)[[1]]) 

#interactive plot
ggplotly(edsvl.plot)

PGLS of log eye diameter vs. log cornea diameter

#PGLS eye diameter vs. cornea diameter
pgls_edcd <- pgls(log10(eye_av) ~ log10(cor_av), 
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edcd)

par(mfrow = c(1, 1))

#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.edcd <- pgls.profile(pgls_edcd, "lambda")
plot(lambda.edcd)

#print model output 
summary(pgls_edcd)

## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(cor_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0069593 -0.0012901 -0.0000271  0.0017229  0.0066368 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.347
##    lower bound : 0.000, p = 6.9782e-05
##    upper bound : 1.000, p = < 2.22e-16
##    95.0% CI   : (0.096, 0.648)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   0.169381   0.012943  13.087 < 2.2e-16 ***
## log10(cor_av) 0.962429   0.021166  45.471 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002278 on 136 degrees of freedom
##   (14 observations deleted due to missingness)
## Multiple R-squared: 0.9383,  Adjusted R-squared: 0.9378 
## F-statistic:  2068 on 1 and 136 DF,  p-value: < 2.2e-16

#plot
edcd.plot <- ggplot(data = av.sal, aes(x = log10(cor_av), y = log10(eye_av), text = Genus_Species)) +
 geom_point(alpha = 0.7) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_edcd)[[2]], intercept = coef(pgls_edcd)[[1]]) 

#interactive plot
ggplotly(edcd.plot)

PGLS of log cornea diameter vs. log cube root of mass

#PGLS cornea diameter vs. cube root of mass
pgls_cdmass <- pgls(log10(cor_av) ~ log10(rootmass_av),
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_cdmass)

par(mfrow = c(1, 1))

#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.cdmass <- pgls.profile(pgls_cdmass, "lambda")
plot(lambda.cdmass)

#print model output 
summary(pgls_cdmass)

## 
## Call:
## pgls(formula = log10(cor_av) ~ log10(rootmass_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.031038 -0.005368 -0.000095  0.004781  0.041013 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.832
##    lower bound : 0.000, p = 8.4377e-15
##    upper bound : 1.000, p = 5.6595e-09
##    95.0% CI   : (0.683, 0.922)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)        0.184679   0.057335  3.2211  0.001566 ** 
## log10(rootmass_av) 0.586098   0.051971 11.2773 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01006 on 150 degrees of freedom
## Multiple R-squared: 0.4588,  Adjusted R-squared: 0.4552 
## F-statistic: 127.2 on 1 and 150 DF,  p-value: < 2.2e-16

#plot
cdmass.plot <- ggplot(data = av.sal, aes(x = log10(rootmass_av), y = log10(cor_av), text = Genus_Species)) +
 geom_point(alpha = 0.7) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_cdmass)[[2]], intercept = coef(pgls_cdmass)[[1]]) 

#interactive plot
ggplotly(cdmass.plot)

PGLS of log cornea diameter vs. log snout-vent length

#PGLS cornea diameter vs. SVL
pgls_cdsvl <- pgls(log10(cor_av) ~ log10(svl_av),
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_cdsvl)

par(mfrow = c(1, 1))

#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.cdsvl <- pgls.profile(pgls_cdsvl, "lambda")
plot(lambda.cdsvl)

#print model output 
summary(pgls_cdsvl)

## 
## Call:
## pgls(formula = log10(cor_av) ~ log10(svl_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.028530 -0.006068  0.001092  0.006296  0.043279 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.872
##    lower bound : 0.000, p = < 2.22e-16
##    upper bound : 1.000, p = 1.4035e-07
##    95.0% CI   : (0.752, 0.942)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   -0.664239   0.120396 -5.5171 1.475e-07 ***
## log10(svl_av)  0.537185   0.051402 10.4507 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01098 on 150 degrees of freedom
## Multiple R-squared: 0.4213,  Adjusted R-squared: 0.4175 
## F-statistic: 109.2 on 1 and 150 DF,  p-value: < 2.2e-16

#plot
cdsvl.plot <- ggplot(data = av.sal, aes(x = log10(svl_av), y = log10(cor_av), text = Genus_Species)) +
 geom_point(alpha = 0.7) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_cdsvl)[[2]], intercept = coef(pgls_cdsvl)[[1]]) 

#interactive plot
ggplotly(cdsvl.plot)

PGLS of log snout-vent length vs. log cube root of mass

#PGLS svl vs. cube root of mass
pgls_svlmass <- pgls(log10(svl_av) ~ log10(rootmass_av),
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

#diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_svlmass)

par(mfrow = c(1, 1))

#Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.svlmass <- pgls.profile(pgls_svlmass, "lambda")
plot(lambda.svlmass)

#print model output 
summary(pgls_svlmass)

## 
## Call:
## pgls(formula = log10(svl_av) ~ log10(rootmass_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0190565 -0.0039221  0.0004095  0.0038784  0.0182303 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.928
##    lower bound : 0.000, p = < 2.22e-16
##    upper bound : 1.000, p = 7.3054e-07
##    95.0% CI   : (0.855, 0.970)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)        1.617161   0.033287  48.582 < 2.2e-16 ***
## log10(rootmass_av) 1.004421   0.025367  39.596 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005716 on 150 degrees of freedom
## Multiple R-squared: 0.9127,  Adjusted R-squared: 0.9121 
## F-statistic:  1568 on 1 and 150 DF,  p-value: < 2.2e-16

#plot
svlmass.plot <- ggplot(data = av.sal, aes(x = log10(rootmass_av), y = log10(svl_av), text = Genus_Species)) +
 geom_point(alpha = 0.7) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(pgls_svlmass)[[2]], intercept = coef(pgls_svlmass)[[1]]) 

#interactive plot
ggplotly(svlmass.plot)

Figures

Here, we plot absolute eye size (cornea diameter) and relative eye size (the residuals of cornea diameter vs. cube root of mass) onto the salamander phylogeny. We will use the ggtree package for plotting. Again, the gilled Pleurodes waltl will not be included in these plots (as it was not included in the PGLS regressions).

First we add PGLS residuals to the dataframe and resort the data to match the species order of the phylogeny.

#Add residuals from PGLS fits to dataframe

#extract pgls residuals 
pglsres.cdmass <- as.data.frame(residuals(pgls_cdmass))

#name residuals
colnames(pglsres.cdmass) <- "pglsres.cdmass" 

#add column for genus_species
pglsres.cdmass$Genus_Species <- rownames(pglsres.cdmass)

#merge residuals with original data by rowname
sal.res <- inner_join(av.sal, pglsres.cdmass, by = "Genus_Species") 

#export dataset for use later
write.csv(sal.res, file = "../Data/Tidy/salamanders_residuals.csv")

#make dataframe for labeling phylogeny figure
sal.phy <- sal.res %>%
  #make column labeling tip labels that match phylogeny
  mutate(tip = Genus_Species) %>%
  #put genus and species in separate columns
  separate(Genus_Species, c("genus", "species"), sep = "_", extra = "drop") %>%
  #add tip labels with regular text
  mutate(labels = as.factor(paste(genus, species, sep = " "))) %>%
  #keep only columns we will use for plotting
  select(tip, labels, cor_av, pglsres.cdmass)

# set row names in dataset to match the tip labels in the tree
row.names(sal.phy) <- sal.phy$tip

#check that phylogeny and data match exactly
name.check(caudatatree, sal.phy)

## [1] "OK"

#resort  dataset to the order of tree tip labels
sal.phy <- sal.phy[caudatatree$tip.label, ] 

#labels
labs <- sal.phy %>% select(tip, labels)

Next, we use ggtree to plot the phylogeny alongside bar plots for absolute cornea diameter and relative conea size.

# Make the phylogeny plot
p <- ggtree(caudatatree) %<+% labs + 
     geom_tiplab(size = 3, aes(label = labels), fontface = 3) + 
   xlim_tree(375) + 
   coord_cartesian(clip = 'off') 

# Make a second plot next to the phylogeny (bars for cornea diameter)
p2 <- facet_plot(p, panel="Cornea diameter (mm)", data=sal.phy, geom=geom_segment, aes(x=0, xend=cor_av, y=y, yend=y), color = "#BDC3BE", size = 3)

# Make a third plot next to the first two (bars for relative cornea size)
p3 <- facet_plot(p2, panel='Cornea ~ mass residuals', data=sal.phy, geom=geom_segment, aes(x=0, xend=pglsres.cdmass, y=y, yend=y), color = "#BDC3BE", size=3) +
  theme_tree2(legend.position=c(.07, .93)) #add scale bars and move legends

#control widths of panels
 gt = ggplot_gtable(ggplot_build(p3))
# gtable_show_layout(gt) # will show you the layout - very handy function
# gt # see plot layout in table format
# gt$layout$l[grep('tree', gt$layout$name)] # you want to find the column specific to panel-2
 gt$widths[7] = 0.5*gt$widths[5] # reduce column 7
 gt$widths[9] = 0.5*gt$widths[5]

#print figure
plot(gt)

#export figure
pdf("../Outputs/tree-figure.pdf", width = 8, height = 22)
grid::grid.draw(gt) # plot with grid draw
dev.off()

## quartz_off_screen 
##                 2